library(sasld)

# -------------------------------------------------------------------
# Sezione 14.1 - L'ANOVA a una via
# -------------------------------------------------------------------
adv <- c(5,1,11,2,8)  # gruppo 1
muz <- c(0,1,4,6,3)   # gruppo 2
cla <- c(13,9,8,15,7) # gruppo 3
y <- c(adv,muz,cla)

(devt <- sum((y - mean(y))^2))

m1 <- mean(adv)
m2 <- mean(muz)
m3 <- mean(cla)
mg <- mean(y)
c(m1,m2,m3,mg)
(devf <- sum((m1-mg)^2+(m2-mg)^2+(m3-mg)^2)*5)

dev1 <- sum((adv - mean(adv))^2)
dev2 <- sum((muz - mean(muz))^2)
dev3 <- sum((cla - mean(cla))^2)
c(dev1,dev2,dev3)
(devr <- dev1 + dev2 + dev3)

(varf <- devf/2)
(varr <- devr/12)

(F <- varf/varr)
pf(F,2,12,lower.tail=FALSE)
qf(0.05,2,12,lower.tail=FALSE)
qf(0.01,2,12,lower.tail=FALSE)


# L'ANOVA a 1 via con la funzione aov
y <- c(adv,muz,cla)
x <- c("A","A","A","A","A",
       "B","B","B","B","B",
       "C","C","C","C","C")

fit <- aov(y ~ x)
anova(fit)


# La verifica degli assunti dell'ANOVA
fit <- aov(y ~ x)
res <- resid(fit)
qqnorm(res)
qqline(res)

levene(y,factor(x))




# -------------------------------------------------------------------
# Sezione 14.2 - Intervalli di confidenza per confrontare coppie di medie
# -------------------------------------------------------------------

# Intervalli di confidenza con il metodo di Fisher
adv <- c(5,1,11,2,8)
muz <- c(0,1,4,6,3)
cla <- c(13,9,8,15,7)
y <- c(adv,muz,cla)
x <- c("A","A","A","A","A",
       "B","B","B","B","B",
       "C","C","C","C","C")
fit <- aov(y ~ x)

m <- tapply(y,x,mean)
n <- table(x)

(MSE <- sum(resid(fit)^2)/fit$df.residual)

#
level <- 0.95
alpha <- 1 - level
amez <- alpha/2
k <- length(m)
kk <- k*(k-1)/2
s <- sqrt(MSE)
glres <- sum(n) - k
out <- matrix(nrow=kk,ncol=10)
kkk <- 0
for (i in 1:(k-1)) {
    for (j in (i+1):k) {
        dlt <- m[i] - m[j]
        tcr <- qt(amez,glres,lower.tail = FALSE)
        me <- tcr*s*sqrt(1/n[i] + 1/n[j])
        lwr <- dlt - me
        upp <- dlt + me
        kkk <- kkk + 1
        out[kkk,1: 6] <- c(i,j,m[i],m[j],n[i],n[j])
        out[kkk,7:10] <- c(dlt,me,lwr,upp)
}
}
colnames(out) <- c("gruppo","gruppo","media","media","n","n",
                   "delta","m.err.","e.inf.","e.sup.")
out


# Esempio del testo
m <- c(8.3, 7.4, 10.4)
n <- c(87, 468,276)
MSE <- 234.8

level <- 0.95
alpha <- 1 - level
amez <- alpha/2
k <- length(m)
kk <- k*(k-1)/2
s <- sqrt(MSE)
glres <- sum(n) - k
out <- matrix(nrow=kk,ncol=10)
kkk <- 0
for (i in 1:(k-1)) {
    for (j in (i+1):k) {
        dlt <- m[i] - m[j]
        tcr <- qt(amez,glres,lower.tail = FALSE)
        me <- tcr*s*sqrt(1/n[i] + 1/n[j])
        lwr <- dlt - me
        upp <- dlt + me
        kkk <- kkk + 1
        out[kkk,1: 6] <- c(i,j,m[i],m[j],n[i],n[j])
        out[kkk,7:10] <- c(dlt,me,lwr,upp)
}
}
colnames(out) <- c("gruppo","gruppo","media","media","n","n",
                   "delta","m.err.","e.inf.","e.sup.")
out


# Intervalli di confidenza con la correzione di Bonferroni
adv <- c(5,1,11,2,8)
muz <- c(0,1,4,6,3)
cla <- c(13,9,8,15,7)
y <- c(adv,muz,cla)
x <- c("A","A","A","A","A",
       "B","B","B","B","B",
       "C","C","C","C","C")
fit <- aov(y ~ x)
m <- tapply(y,x,mean)
n <- table(x)
(MSE <- sum(resid(fit)^2)/fit$df.residual)
#
level <- 0.95
alpha <- 1 - level
amez <- alpha/2
k <- length(m)
kk <- k*(k-1)/2
s <- sqrt(MSE)
glres <- sum(n) - k
out <- matrix(nrow=kk,ncol=10)
kkk <- 0
for (i in 1:(k-1)) {
    for (j in (i+1):k) {
        dlt <- m[i] - m[j]
        tcr <- qt(amez/kk,glres,lower.tail = FALSE) # <--------------
        me <- tcr*s*sqrt(1/n[i] + 1/n[j])
        lwr <- dlt - me
        upp <- dlt + me
        kkk <- kkk + 1
        out[kkk,1: 6] <- c(i,j,m[i],m[j],n[i],n[j])
        out[kkk,7:10] <- c(dlt,me,lwr,upp)
}
}
colnames(out) <- c("gruppo","gruppo","media","media","n","n",
                   "delta","m.err.","e.inf.","e.sup.")
out


# Intervalli di confidenza con il metodo di Tukey
fit <- aov(y ~ x)
TukeyHSD(fit)




# -------------------------------------------------------------------
# Sezione 14.3 - ANOVA e regressione
# -------------------------------------------------------------------
adv <- c(5,1,11,2,8); muz <- c(0,1,4,6,3); cla <- c(13,9,8,15,7)
y <- c(adv,muz,cla)
x <- c("A","A","A","A","A","B","B","B","B","B","C","C","C","C","C")
x <- factor(x)
contrasts(x)

x <- relevel(x, ref="C")
contrasts(x)

fit <- lm(y ~ x)
summary(fit)

confint(fit)

# per ottenere l'i.c. mancante (A-B)
x <- relevel(x, ref="B")
fit <- lm(y ~ x)
confint(fit)




# -------------------------------------------------------------------
# Sezione 14.4 - ANOVA a due vie
# -------------------------------------------------------------------
y <- c(137,158,139,166,155,
       164,125,141,144,122,
       150,151,120,157,122,
       124,106,137, 87,109)/10

fertilizer <- c("H","H","H","H","H","H","H","H","H","H",
                "L","L","L","L","L","L","L","L","L","L")
manure     <- c("H","H","H","H","H","L","L","L","L","L",
                "H","H","H","H","H","L","L","L","L","L")

# ANOVA a due vie (solo effetti principali)
fit <- aov(y ~ fertilizer+manure)
anova(fit)

# La regressione con le variabili indicatrici in una ANOVA a due vie
fertilizer <- factor(fertilizer,levels=c("L","H"))
manure <- factor(manure,levels=c("L","H"))
fit <- lm(y ~ fertilizer+manure)
summary(fit)

fv <- fitted(fit)
tapply(fv,list(manure,fertilizer),mean)
tapply(y,list(manure,fertilizer),mean)


# per avere lo stesso output del pdf
fertilizer <- factor(fertilizer,levels=c("H","L"))
manure <- factor(manure,levels=c("H","L"))
#
tapply(y,fertilizer,mean)
tapply(fv,fertilizer,mean)

tapply(y,manure,mean)
tapply(fv,manure,mean)


# ridefinisco i livelli per avere lo stesso output del pdf
fertilizer <- factor(fertilizer,levels=c("H","L"))
manure <- factor(manure,levels=c("H","L"))
fit <- lm(y ~ fertilizer+manure)
confint(fit)


# ANOVA a due vie (effetti principali e interazione)
fit <- aov(y ~ fertilizer+manure+fertilizer:manure)
anova(fit)

fit <- lm(y ~ fertilizer+manure+fertilizer:manure)
summary(fit)

fv <- fitted(fit)
tapply(fv,list(manure,fertilizer),mean)


# ANOVA a due vie in un disegno sbilanciato
# considero solo i dati dell'esempio nel testo
data(gss)

nok <- which(gss$race == "other"); gss <- gss[-nok,]
gss$race <- as.character(gss$race)
y <- as.numeric(gss$polviews) - 1; gss$y <- y
nok <- which(is.na(y)); gss <- gss[-nok,]
#
table(gss$sex,gss$race)

#
fit <- aov(y ~ sex+race, data=gss)
anova(fit)
#
fit <- aov(y ~ race+sex, data=gss)
anova(fit)

fit <- aov(y ~ sex*race, data=gss)
anova(fit)
#
fit <- aov(y ~ race*sex, data=gss)
anova(fit)

fit <- aov(y ~ sex+race, data=gss)
drop1(fit, test="F")